Airy function was firstly proved to be an exact solution of the Schrödinger equation by Berry and Balazs et al. in 1979.^[1] But ideal Airy beam carries infinite energy, it is not physically realizable. Since Siviloglou first proved the existence of the finite energy Airy beam with truncation parameter in the experiment in 2007,^[2,3] the Airy beam has attracted wide attention because of its nondiffraction, self-accelerating and self-healing properties,^[4–6] which come into being in various potential applications, including optical particle control,^[7] electron acceleration,^[8,9] the generation of light bullets and plasma channel,^[10,11] and so on.

Recently, the research of spatiotemporal optical wave packets or light bullets has also been a fascinating field, especially the application of the Airy function in optical wave packet or light bullet generation:^[11–17] the nondispersion and nondiffraction Airy–Bessel wave packets which combine one Airy function in the time domain and Bessel function in space domain have been reported as versatile linear light bullets by Chong et al.;^[11] the generation mechanism of Airy–Bessel wave packets in free space have been researched by Ren et al.;^[12] the realization of spatiotemporal Airy light bullets has been demonstrated by Abdollahpour et al. through combining a spatial Airy function and a temporal Airy function.^[13] In addition, in the help of the Airy function, the (3+1)-dimensional ((3+1)D) localized Airy elegant Laguerre–Gaussian and Airy elegant Hermite–Gaussian wave packets have been studied respectively by Wu et at.^[15] and Deng et al.^[17]

On the other hand, to explore the properties and applications of paraxial optical wave packets with elliptic geometry has also been an interesting study field. The elegant Ince–Gaussian (eIG) function that can be regarded as the arithmetic product of Ince polynomials with complex arguments and a Gaussian function has been theoretically demonstrated to be the solution of the paraxial wave equation in elliptic coordinates,^[18] and is also complementary to the Ince–Gaussian functions that are the third complete family of exact solutions of the paraxial wave equation in free space.^[19–22] So, the eIG function will be a good choice for constructing wave packets with elliptic geometry. In addition, the Airy helical elegant Ince Gaussian (AiHeIG) wave packets with orbital angular momentum have not been studied yet. Therefore, it will also be interesting to study what will happen when Airy elegant Ince–Gaussian (AieIG) wave packets carry orbital angular momentum.

In this paper, we study a class of (3+1)D AieIG and AiHeIG wave packets in free space obtained with the help of the modulated Airy pulse and Ince polynomials with complex arguments. We demonstrate the numerical experiments of these (3+1)D wave packets at Z = 0, T = 0 and analyze the evolution properties and the energy conversion process of these self-accelerating spatiotemporal wave packets. The results show that these wave packets will self-decelerate along the time domain and their shapes and energy distributions will change during propagation. Lastly, we analyze the radiation forces and 3-dimensional (3D) intensity and phase of these wave packets, which will help us understand the structure of the kind of elliptic wave packets more comprehensively.

The Ince polynomials with complex arguments can be used to describe (3+1)D finite energy wave packets in the presence of diffraction and dispersion when combined with other nondiffracting field configurations, for example Airy pulse. In this case, the whole wave packets in the linear spatiotemporal domain within the framework of the paraxial approximation evolution obey the equation:^[3,23] 1 where Ψ(X,Y,Z,T) denotes the slowly varying complex envelope of the (3+1)-dimensional optical field; X, Y are the dimensionless transverse coordinates that are scaled by the waist spot w₀; T is the dimensionless time coordinate that is scaled by the temporal scaling parameter t₀; Z = z/Z_R represents the normalized evolution distance, and is the Rayleigh length, λ₀ is the vacuum wavelength.

In general, the scientific workers always expect the usual Gaussian factor exp[(ik/2q)r²] to be an essential part of the eigenfunction when solving the paraxial wave equation, that is the reason why we assume our proposed AieIG wave packets as the form: 2 where A, E, N are complex functions and Ψ_G = exp[(ik/2q)r²] is the Gaussian factor; r is the radius; c = c(Z) = 1/2i(Z − i) is a complex parameter. In the transverse plane perpendicular to Z, the complex elliptic coordinate transformation is defined as: X = f(Z) cosh ξ cos η, Y = f(Z) sinh ξ sin η, and Z = Z, where ξ and η are respectively the radial and angular elliptic variables, f(Z) = f₀/(w₀c(Z)^1/2) and f₀ is the semi-focal separation at the waist plane Z = 0, which is the only real elliptical coordinate plane in the whole spatial domain because c(0) ≡ 1/2. Putting Eq. (2) into Eq. (1), one can obtain the separated three equations: 3 4 5 where a and p are the separation constants; and ε is the ellipticity parameter, which is decided by the waist spot w₀ and the semi-focal separation f₀, exactly, . The ellipticity ε, waist spot w₀ and semi-focal separation f₀ are all the physically important parameters for describing the transverse structure of the AieIG wave packets. The ellipticity parameter ε adjusts the ellipticity of modes, while the parameters w₀ and f₀ scale their physical size.^[20,21]

For finding the solution of Eq. (3), we let A(Z,T) = (−i/Z − i)^p+1A₁(Z,T), then equation (3) can be transformed into 6 Equation (6) is also known as the normalized (1+1)D paraxial diffraction equation, which can be solved by the Fourier-transform method and the self-accelerating Airy beam solution of Eq. (6) has been discussed extensively.^[1–3] So, we can obtain easily the solution of Eq. (6) that can be written as 7 where Ai(·) is the Airy function, and 0 < α < 1 is the decay parameter. Equation (7) is in keeping with the expression for spatial Airy beam^[3] if we use the spatial coordinate s instead of the temporal coordinate T.

Then we obtain the solution of Eq. (3) as follows: 8

Equation (8) shows the intensity profile of the Airy pulse modulated by the transmission distance Z and the separation parameter p, whose peak intensity decreases with T. Figure 1(a) displays the normalized intensity profile at several different evolution distances. As seen in Fig. 1(a), the normalized peak intensity of the modulated Airy pulse decreases fleetingly with the increasing of transmission distance Z. The propagation dynamic of the pulse is depicted in Fig. 1(b), where we can see the propagation process more clearly, including the self-accelerating process and the propagation trajectory of the pulse.

Fig. 1.

	Figure Option ViewDownloadNew Window
	Fig. 1. (color online) (a) The normalized intensity profile of the modulated Airy pulse in different evolution distances, (b) numerical simulation of the propagation dynamic of the modulated Airy pulse as a function of the normalized distance. The parameters are α = 0.1, p = 4.

Equation (4) is a special case of the most ordinary Hill equation, also known as the Ince equation.^[22,24–26] One can derive Eq. (5) from Eq. (4) by replacing ξ with −iη, and vice versa. The interconversion is advantageous because it makes one get the radial solutions E(ξ) from angular solutions N(η) more easily. On the other hand, solutions of Eq. (5) are divided into the two cases that are known as the even and odd Ince polynomials with order p and degree m, generally denoted as and respectively, where p ≥ m, and the parameters (p,m) must have the same parity, i.e, (−1)^p−m = 1 must be obeyed. The minimum m is 0 for even cases, and 1 for odd cases. In the process of searching solutions, only when the same parities in both ξ and η satisfy continuity in the whole space, we can obtain the right solution. The general expressions of the even and odd eIG functions can be written respectively as: 9 10 where C_e and S_e are normalization constants, the superscripts e and o denote respectively even and odd cases.

Substituting the Eqs. (8)–(10) into the Eq. (2), we can obtain the solution of the Eq. (1) named the self-accelerating even or odd AieIG wave packets, which are expressed respectively as: 11 12

The self-accelerating AiHeIG wave packets can also be regarded as a linear combination of the self-accelerating even and odd AieIG wave packets. Here, we give the common expression of the self-accelerating AiHeIG wave packets as follows: 13 where the sign H⁺ denotes the positive helixes. Note that m > 0, because odd Ince polynomials do not include the case m = 0.

Figure 2 shows the normalized transverse intensity and phase distributions of some low-order AiHeIG wave packets at z = 0, T = 0. As shown in Fig. 2, these input wave packets can be described by elliptic rings and vortices that are physically separated in-line, which are controlled by the parameters p and m, namely elliptic rings N = 1 + (p − m)/2, and the number of vortices is m. The initial input wave packets with the same m have the similar structure. On the other hand, the parameter p determines the cylinder number of wave packets. More exactly, the parameter p adds 2 and the cylinder number will add 1 in terms of wave packets with the same m [see Fig. 2]. It is also noted that the intensity of wave packets' lateral cylinder is smaller than that in the inner at T = 0, Z = 0. Usually, we just can see three circles at most, but the phase distribution can tell us accurately how many circles they have.

Fig. 2.

	Figure Option ViewDownloadNew Window
	Fig. 2. (color online) Normalized transverse intensity distributions of some AiHeIG wave packets, ε = 2, Z = 0, and T = 0. Plots in the bottom row of every mode correspond to the normalized phase structures.

We have obtained the accurate analytical solutions of the (3+1)D self-accelerating spatiotemporal wave packets in free space from (3+1)D potential-free Schrödinger equation, namely, Eqs. (11)–(13). Evidently, the varieties of the even, odd and helical AieIG wave packets can be constructed by adjusting the indices (p,m). Therefore, we divide these wave packets into two typical cases that p ≠ m and p = m and investigate their properties respectively. In these two cases, we study the different situations of the wave packets with varying parameters and discover many exciting phenomena. In addition, we also give the discussion of the 3D intensity and phase of these wave packets.

3.1. Case p ≠ m

Firstly, we discuss the case p ≠ m. For simplicity, we just discuss the case ε = 2 here and take p = 4, m = 2 as an example. Figures 3(a1)–3(c1) are respectively the normalized transverse intensity patterns of even, odd, and helical AieIG wave packets at Z = 0, T = 0. As we expected, the initial input AiHeIG_4,2 wave packets can be described by two elliptical rings and two coaxial vortexes, but the lateral rings is too faint to be observed [see Fig. 3(c1)]. Figures 3(a2)–3(c2) are the corresponding phase distributions of the initial input wave packets. The normalized interference intensity distributions are obtained by calculating the off-axis interference patterns between the complex amplitude profiles of even, odd and helical AieIG wave packets at the Z = 0 plane and a plane wave when T = 0 [see Figs.3(a3)–3(c3)]. For the numerical experimental generation, initial wave packets have been launched to reconstruct the off-axis computer-generated holograms [see Figs. 3(a4)–3(c4)] of the desired wave packets.^[27] Figures 3(a5)–3(c5) numerically simulate the normalized transverse intensity distributions of the even, odd and helical AieIG wave packets at Z = 1, T = 0, as contrast. The normalized transverse intensity patterns at Z = 1 obtained by numerical experiment are shown as Figs. 3(a6)–3(c6). As we can see, our results of accurately numerical experiment are in good agreement with our analytic results. In addition, the intensity distributions vary with the evolution distance. Due to the influence of vortices, we can see the rotation effect in the process of AiHeIG wave packets' evolution. It seems that there is a centrifugal force that pulls these AiHeIG wave packets to produce this phenomenon.

The corresponding isosurface intensity plots of these (3+1)D self-accelerating AieIG wave packets are shown in Fig. 4. They are accelerated in the temporal domain. Due to the effect of the modulated Airy pulse, the profiles of these wave packets are not distinguished in T axis with the increasing of the evolution distance. That is to say, the pulse stack occurs in the process of these (3+1)D wave packets' evolution. In addition, for even AieIG wave packets, the center energy decreases rapidly and transfers to sidelobes in the evolution process of these wave packets so that we can see more lateral sidelobes as shown in Figs. 4(a1)–4(a3). A result that can be observed is that the pulse stack of the wave packets' center is decreasing and that of marginal wave packets is increasing with the increasing of the evolution distance. However, for odd AieIG_4,2 wave packets, due to the energy balance of these wave packets' every part at initial input plane, there is no obvious energy conversion in every part [see Figs. 3(b1), 3(b5), 3(b6), and Figs. 4(b1)–4(b3)], which is obviously different from the even AieIG_4,2 wave packets. The most obvious feature of the odd AieIG_4,2 is that each lobe of the wave packets is equably extending along the x axis and y axis with the increasing of the evolution distance. For AiHeIG wave packets, we still see the obvious rotation phenomenon which will also be accompanied by the energy conversion.

Fig. 3.

Figure Option T = 0. (a1)–(a6) the even AieIG wave packets; (b1)–(b6) the odd AieIG wave packets; (c1)–(c6) the helical AieIG wave packets. (a1)–(c1) numerical simulation of the normalized intensity distributions at Z = 0; (a2)–(c2) numerical simulation of the normalized phase distributions at Z = 0; (a3)–(c3) the normalized interference intensity of the initial generated beams and a plane wave; (a4)–(c4) the computer-generated holograms; (a5)–(c5) numerical simulation of the normalized transverse intensity distributions at Z = 1, and (a6)–(c6) the corresponding numerical experimental record. The other parameters are chosen as p = 4, m = 2, and ε = 2. '>ViewDownloadNew Window

T = 0. (a1)–(a6) the even AieIG wave packets; (b1)–(b6) the odd AieIG wave packets; (c1)–(c6) the helical AieIG wave packets. (a1)–(c1) numerical simulation of the normalized intensity distributions at Z = 0; (a2)–(c2) numerical simulation of the normalized phase distributions at Z = 0; (a3)–(c3) the normalized interference intensity of the initial generated beams and a plane wave; (a4)–(c4) the computer-generated holograms; (a5)–(c5) numerical simulation of the normalized transverse intensity distributions at Z = 1, and (a6)–(c6) the corresponding numerical experimental record. The other parameters are chosen as p = 4, m = 2, and ε = 2. '>

Fig. 3. (color online) Numerical demonstrations of the even, odd, and helical AieIG wave packets' evolution with different distances when T = 0. (a1)–(a6) the even AieIG wave packets; (b1)–(b6) the odd AieIG wave packets; (c1)–(c6) the helical AieIG wave packets. (a1)–(c1) numerical simulation of the normalized intensity distributions at Z = 0; (a2)–(c2) numerical simulation of the normalized phase distributions at Z = 0; (a3)–(c3) the normalized interference intensity of the initial generated beams and a plane wave; (a4)–(c4) the computer-generated holograms; (a5)–(c5) numerical simulation of the normalized transverse intensity distributions at Z = 1, and (a6)–(c6) the corresponding numerical experimental record. The other parameters are chosen as p = 4, m = 2, and ε = 2.

Fig. 4.

	Figure Option ViewDownloadNew Window
	Fig. 4. (color online) Isosurface intensity plots of the self-accelerating even, odd, and helical AieIG wave packets with p = 4, m = 2, and ε = 2. (a1)–(a3) the even case; (b1)–(b3) the odd case; (c1)–(c3) the helical case; (a1)–(c1) Z = 0; (a2)–(c2) Z = 1; (a3)–(c3) Z = 2.

In this part, we will discuss the case p = m and take p = 4, m = 4 as an example. When the ellipticity ε takes the limit, the wave packets’ pattern changes, so the number of coaxial vortexes changes as well.

As shown in the following Fig. 5, when ε → 0, these separated vortexes spatially collapse to a single vortex with the total charge of m and the AiHeIG wave packets will transit into the Airy helical Laguerre–Gaussian wave packets; when the ellipticity ε = 2 or ε → ∞, the coaxial vortexes turn up. However, when ε = 2, namely AiHeIG wave packets, the number of the vortexes is decided by the parameter m; when ε → ∞, namely Airy helical elegant Hermite–Gaussian wave packets, the vortex number is decided by the number of the optical singular points that the real and imaginary parts of wave packets are zero simultaneously.

From Figs. 5(a3)–5(c3) and 5(a4)–5(c4), we can see that our results of accurately numerical experiment are in good agreement with our analytic results. The clearly temporal evolution processes are shown in Figs. 5(a5)–5(c5). We still can see the rotation phenomenon by comparing with Figs. 5(a1)–5(c1), respectively.

Fig. 5.

Figure Option ViewDownloadNew Window

Fig. 5. (color online) Numerical demonstrations of the AiHeIG wave packets with different ellipticity parameters, distances Z and time T: the top row with the case ε → 0; the middle row with the case ε = 2; the bottom row with the case ε → ∞. (a1)–(c1) With the intensity distribution at Z = 0, T = 0; (a2)–(c2) with the corresponding phase distribution at Z = 0, T = 0; (a3)–(c3) the corresponding intensity distribution at Z = 1, T = 0; (a4)–(c4) the corresponding numerical experimental record intensity distribution at Z = 1, T = 0; (a5)–(c5) snapshots describing the evolution of the self-accelerating AieIG wave packets at Z = 1.

3.2. The three-dimensional intensity and phase of AiHeIG wave packets

In order to study the characteristic of the AieIG wave packets in free space better, we also give the three-dimensional intensity and phase distributions of several helical AieIG wave packets.

Figures 6(a1)–6(a3) show respectively the three-dimensional propagation of AiHeIG_4,2, AiHeIG_3,3, and AiHeIG_4,4 wave packets which cut off the section Y greater than 0. From Figs. 6(a1)–6(a3), we can see more clearly the process that the central energy of these wave packets is gradually transferred to nearby sidelobes until disappeared and the evolution process of vortexes. In addition, as we can see, the number of vortexes increases with the parameter m. All the wave packets trend to two different route propagations finally. Figures 6(b1)–6(b3) show the corresponding normalized phase distribution on section, which is important to restructure these wave packets and can also help us understand the characteristics of these wave packets well.

Fig. 6.

	Figure Option ViewDownloadNew Window
	Fig. 6. (color online) The isosurface and section plot of the AiHeIG wave packets during propagation. (a) the corresponding normalized intensity distributions; (b) the corresponding normalized phase distributions. The left column with the case (4,2); the middle column with the case (3,3); the right column with the case (4,4).

The radiation forces play an important role in analysis optical phenomena. The gradient force and the scattering force are considered as two kinds of common radiation forces. Assume that there is a particle in the space whose radius r₀ is sufficiently small compared with the wavelength, and can be seen as a point dipole in the light field. With the help of the vector identity and the solution of the Maxwell equations, the gradient force on a Rayleigh dielectric particle^[28] can be written as 14 where is the electric dipole moment of the particle,^[29,30] n is the refractive index of the particulate, ε₀ is the permittivity of vacuum. Assuming the particle is in stable state, the gradient force will be the time average,^[28] namely: 15 where c is the light velocity,

The change of the electromagnetic momentum will influence the scattering of light, with which the scattering force is associated. The scattering force^[28] can be expressed as 16 where C_pr0 is the radiation pressure section of the particle, e_z is a unity vector along the beam propagation. Considering the isotropy of the particle, the radiation pressure section of the particle is equal to the scattering force section of the scatterer, namely: 17

Figure 7 shows the numerically simulated transverse gradient forces of several AieIG wave packets. The particle parameters are chosen as n = 1.2, r₀ = 60 nm. Compared with Figs. 3(a1), 3(a5), and 3(a6), we can find the relation of the amplitude distribution of gradient forces with that of the optical intensity. Generally, the gradient forces will be 0 at the position where optical intensity has its peak, so we can see some vortexes in those positions. As shown in Figs. 7(c1)–7(c3), the amplitude distributions of gradient forces have the nested double light intensity structure and they also rotate with wave packets' evolution. Figures 7(a4)–7(c4) show the corresponding gradient forces along with x axis. The sign of the forces is determined by the direction of the radiation force. The peak of gradient forces of the even AieIG_3,1 and AieIG_4,2 diminishes fast with the increasing of distance Z, nevertheless, that of the helical wave packets reduces slower.

Fig. 7.

Figure Option ViewDownloadNew Window

Fig. 7. (color online) The gradient force of several self-accelerating AieIG wave packets at T = 0. (a1)–(a3) with the even (3,1) case; (b1)-(b3) with the even (4,2) case; (c1)–(c3) with the helical (4,4) case; (d1)–(d3) with the normalized peak of the gradient force. (a1)–(c1) with the case Z = 0; (a2)–(c2) with the case Z = 1; (a3)–(c3) with the case Z = 2; (a4)–(d4) the gradient force of even (3,1), even (4,2) and helical (4,4) along with x-axis.

The scattering force is in proportion to light intensity in terms of constant refractive index particulate. Figure 8 shows the scattering force of several AieIG wave packets propagating at the same plane as those in Fig. 7. It is not difficult to find that there are many similar behaviors between the scattering force and the gradient force when comparing with Fig. 8. The maximum of the scattering force of the odd AieIG_3,1 wave packets extends along longitudinal direction but that of the odd AieIG_3,3 wave packets extends along transverse direction and transmits into other lobes. The normalized peak diminishes and keeps the symmetrical structure with the increase of the distance Z. In addition, due to the proportional relation of scattering force and light intensity, we can also find the energy of the heart of these wave packets decreases and shifts to other positions.

Fig. 8.

Figure Option ViewDownloadNew Window

Fig. 8. (color online) The scattering force of several self-accelerating AieIG wave packets at T = 0. (a1)–(a3) with the odd (3,1) case; (b1)–(b3) with the odd (3,3) case; (c1)–(c3) with the helical (3,3) case; (d1)–(d3) with the normalized peak of the scattering force. (a1)–(d1) with the case Z = 0; (a2)–(d2) with the case Z = 1; (a3)–(d3) with the case Z = 2. All the color scales have been amplified 1021.

In conclusion, we have constructed the even, odd and helical AieIG wave packets' solutions from the (3+1)D potential-free Schrödinger equation with the help of the Airy function and Ince polynomials, and demonstrated the numerical experiment of these AieIG wave packets. Our results show different (3+1)D spatiotemporal wave packets can be presented by selecting different values of the three relevant parameters: indices (p,m) and the ellipticity ε. The controlled behaviors of these (3+1)D spatiotemporal wave packets have been displayed in this paper. A conclusion can be drawn that the different photoelastics, pulse stack, boundary, elliptical rings, physically separated in-line vortices, energy conversion, and wave packets rotation can be controlled by adjusting the ellipticity ε, the evolution distance Z and the mode numbers (p,m). In addition, we also find a very interesting phenomenon that if the energy distribution between every part of these wave packets is uneven at the input plane, the energy will be transferred between every part in the process of transmission. Lastly, we have also analyzed the radiation forces and given 3D intensity and phase distributions that can help us understand the structure of the kind of elliptic wave packets more comprehensively.